Confounder age (Igf & animal welfare)

Author

Anja Eggert & Anne-Marie Galow

Published

March 3, 2025

R Libraries

library(tidyverse)     # tidy universe
library(kableExtra)    # html-table
library(patchwork)     # combine plots
library(lmerTest)      # mixed model
library(car)           # ANOVA
library(performance)   # model performance
set.seed(1989)

Data

Read data

Data processing done in file “igf-biomarker-summary-stats.qmd”.

load("./data/data-processed.RData")

Statistical model: general design

With a simple model structure we solely test for a confounder effect of insemination age. We include a random animal effect and also the nested structure. Each sow gave birth 1, 2 or 3 times, so these events are nested within sows. We fit a linear mixed model with the lmerTest package in R.

contr = lmerControl(optimizer   = "bobyqa",
                    optCtrl     = list(maxfun = 10000000),
                    calc.derivs = FALSE)

Cortisol serum

Model

mod.cort1 <- lmerTest::lmer(log(value) ~ 
                              insem.age +
                              # random intercept for sows and for each litter within sow
                              (1 | sow/litter.no),
                            data    = dat.l |> 
                              dplyr::filter(parameter == "cort.ser") |> 
                              droplevels() |>
                              drop_na(value),
                            REML    = TRUE,
                            control = contr)
summary(mod.cort1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "cort.ser")),  
    value)
Control: contr

REML criterion at convergence: 25.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.88991 -0.48347 -0.01562  0.55338  1.54576 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.01358  0.1166  
 sow           (Intercept) 0.02233  0.1494  
 Residual                  0.05136  0.2266  
Number of obs: 39, groups:  litter.no:sow, 20; sow, 14

Fixed effects:
              Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)  3.7595871  0.2312843  8.5866965  16.255 9.53e-08 ***
insem.age   -0.0011780  0.0006513  8.0024430  -1.809    0.108    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.965
round(drop1(mod.cort1, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age  0.168   0.168     1 8.002   3.271  0.108
car::Anova(mod.cort1,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
           Chisq Df Pr(>Chisq)  
insem.age 3.2711  1    0.07051 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.cort1,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
               F Df Df.res Pr(>F)
insem.age 2.7669  1 9.0492 0.1304

Model diagnostics

performance::check_model(mod.cort1)

Plot

plot.cort1 <- dat.l  |> 
  dplyr::filter(parameter == "cort.ser") |> 
  droplevels() |>
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 60), 
                     breaks = seq(0, 60, 20),
                     minor_breaks = seq(0, 60, by = 2) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "Cortisol (serum) [ng/ml]",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.cort1

Cortisol saliva

Model

mod.cort2 <- lmerTest::lmer(log(value) ~ 
                              insem.age +
                              # random intercept for sows and for each litter within sow
                              (1 | sow/litter.no),
                            data    = dat.l |> 
                              dplyr::filter(parameter == "cort.sal") |> 
                              droplevels() |>
                              drop_na(value),
                            REML    = TRUE,
                            control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.cort2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "cort.sal")),  
    value)
Control: contr

REML criterion at convergence: 115.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.18019 -0.60262  0.02342  0.75178  1.69857 

Random effects:
 Groups        Name        Variance  Std.Dev. 
 litter.no:sow (Intercept) 5.079e-17 7.127e-09
 sow           (Intercept) 0.000e+00 0.000e+00
 Residual                  8.584e-01 9.265e-01
Number of obs: 39, groups:  litter.no:sow, 20; sow, 14

Fixed effects:
              Estimate Std. Error         df t value Pr(>|t|)  
(Intercept)  1.8128656  0.6986073 37.0000000   2.595   0.0135 *
insem.age   -0.0004269  0.0019688 37.0000000  -0.217   0.8295  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.977
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.cort2, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age   0.04    0.04     1    37   0.047   0.83
car::Anova(mod.cort2,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
          Chisq Df Pr(>Chisq)
insem.age 0.047  1     0.8284
car::Anova(mod.cort2,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
               F Df Df.res Pr(>F)
insem.age 0.0402  1 11.608 0.8446

Model diagnostics

performance::check_model(mod.cort2)

Plot

plot.cort2 <- dat.l  |> 
  dplyr::filter(parameter == "cort.sal") |> 
  droplevels() |>
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 30), 
                     breaks = seq(0, 30, 10),
                     minor_breaks = seq(0, 30, by = 2) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "Cortisol (saliva) [ng/ml]",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.cort2

IGF bioactivity

Model

mod.bioact <- lmerTest::lmer(log(value) ~ 
                               insem.age +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "bioact.ser") |> 
                               drop_na(value) |> 
                               dplyr::filter(time != "30dpc") |> 
                               droplevels(),
                            REML    = TRUE,
                            control = contr)
summary(mod.bioact)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
   Data: droplevels(dplyr::filter(drop_na(dplyr::filter(dat.l, parameter ==  
    "bioact.ser"), value), time != "30dpc"))
Control: contr

REML criterion at convergence: 157.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.41286 -0.33591  0.05923  0.42814  2.23561 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.2888   0.5374  
 sow           (Intercept) 0.1375   0.3708  
 Residual                  0.1364   0.3694  
Number of obs: 80, groups:  litter.no:sow, 40; sow, 15

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)  5.391861   0.337537 33.159135  15.974   <2e-16 ***
insem.age   -0.001524   0.000716 26.031980  -2.128   0.0429 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.915
round(drop1(mod.bioact, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)  
insem.age  0.618   0.618     1 26.032   4.529  0.043 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.bioact,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
           Chisq Df Pr(>Chisq)  
insem.age 4.5294  1    0.03332 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.bioact,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
               F Df Df.res  Pr(>F)  
insem.age 4.4425  1  27.39 0.04434 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.bioact)

Plot

plot.bioact <- dat.l  |> 
  dplyr::filter(parameter == "bioact.ser") |> 
  drop_na(value) |> 
  dplyr::filter(time != "30dpc") |> 
  droplevels() |>
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 600), 
                     breaks = seq(0, 600, 200),
                     minor_breaks = seq(0, 600, by = 50) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "IGF bioactivity [ng/ml]",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.bioact

IGF1 (serum)

Model

mod.igf1 <- lmerTest::lmer(log(value) ~ 
                              insem.age +
                              # random intercept for sows and for each litter within sow
                              (1 | sow/litter.no),
                            data    = dat.l |> 
                              dplyr::filter(parameter == "igf1.ser") |> 
                              droplevels() |>
                              drop_na(value),
                            REML    = TRUE,
                            control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.igf1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igf1.ser")),  
    value)
Control: contr

REML criterion at convergence: 102.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6324 -0.6741  0.1403  0.6757  1.5000 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.00441  0.0664  
 sow           (Intercept) 0.00000  0.0000  
 Residual                  0.45051  0.6712  
Number of obs: 44, groups:  litter.no:sow, 22; sow, 15

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)  3.639462   0.471424 20.000005   7.720 2.01e-07 ***
insem.age    0.003985   0.001361 20.000005   2.928  0.00831 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.976
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.igf1, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)   
insem.age  3.863   3.863     1    20   8.574  0.008 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.igf1,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
           Chisq Df Pr(>Chisq)   
insem.age 8.5744  1   0.003409 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.igf1,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
               F Df Df.res  Pr(>F)  
insem.age 7.4073  1 13.461 0.01702 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.igf1)

Plot

plot.igf1 <- dat.l  |> 
  dplyr::filter(parameter == "igf1.ser") |> 
  droplevels() |>
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 500), 
                     breaks = seq(0, 500, 100),
                     minor_breaks = seq(0, 500, by = 20) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "IGF1 (serum) [ng/ml]",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.igf1

IGF2 (serum)

Model

mod.igf2 <- lmerTest::lmer(log(value) ~ 
                              insem.age +
                              # random intercept for sows and for each litter within sow
                              (1 | sow/litter.no),
                            data    = dat.l |> 
                              dplyr::filter(parameter == "igf2.ser") |> 
                              droplevels() |> 
                              drop_na(value),
                            REML    = TRUE,
                            control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.igf2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igf2.ser")),  
    value)
Control: contr

REML criterion at convergence: 61.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2131 -0.8285  0.2732  0.6983  1.3886 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.0000   0.0000  
 sow           (Intercept) 0.0000   0.0000  
 Residual                  0.1937   0.4402  
Number of obs: 40, groups:  litter.no:sow, 20; sow, 14

Fixed effects:
              Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)  4.451e+00  3.230e-01  3.800e+01  13.778 2.33e-16 ***
insem.age   -7.232e-05  9.342e-04  3.800e+01  -0.077    0.939    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.977
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.igf2, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age  0.001   0.001     1    38   0.006  0.939
car::Anova(mod.igf2,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
          Chisq Df Pr(>Chisq)
insem.age 0.006  1     0.9383
car::Anova(mod.igf2,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
              F Df Df.res Pr(>F)
insem.age 0.005  1   12.7 0.9448

Model diagnostics

performance::check_model(mod.igf2)

Plot

plot.igf2 <- dat.l  |> 
  dplyr::filter(parameter == "igf2.ser") |> 
  droplevels() |> 
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 200), 
                     breaks = seq(0, 200, 50),
                     minor_breaks = seq(0, 200, by = 10) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "IGF2 (serum) [ng/ml]",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.igf2

IGFBP2 (serum)

Model

mod.igfbp2 <- lmerTest::lmer(log(value) ~ 
                               insem.age +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "igfbp2.ser") |> 
                               droplevels() |> 
                               drop_na(value),
                             REML    = TRUE,
                             control = contr)
summary(mod.igfbp2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp2.ser")),  
    value)
Control: contr

REML criterion at convergence: 72.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.89723 -0.49844  0.00183  0.48081  2.04124 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.06597  0.2568  
 sow           (Intercept) 0.07302  0.2702  
 Residual                  0.05823  0.2413  
Number of obs: 72, groups:  litter.no:sow, 36; sow, 15

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept) 5.946e+00  1.871e-01 3.105e+01  31.773   <2e-16 ***
insem.age   8.381e-04  4.006e-04 2.365e+01   2.092   0.0473 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.881
round(drop1(mod.igfbp2, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)  
insem.age  0.255   0.255     1 23.649   4.378  0.047 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.igfbp2,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
           Chisq Df Pr(>Chisq)  
insem.age 4.3778  1    0.03641 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.igfbp2,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
               F Df Df.res Pr(>F)  
insem.age 4.2439  1 24.209 0.0503 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.igfbp2)

Plot

plot.igfbp2 <- dat.l  |> 
  dplyr::filter(parameter == "igfbp2.ser") |> 
  droplevels() |> 
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 1250), 
                     breaks = seq(0, 1250, 250),
                     minor_breaks = seq(0, 1250, by = 100) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "IGFBP2 (serum) [ng/ml]",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.igfbp2

IGFBP3 (serum)

Model

mod.igfbp3 <- lmerTest::lmer(log(value) ~ 
                               insem.age +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "igfbp3.ser") |> 
                               droplevels() |> 
                               drop_na(value),
                             REML    = TRUE,
                             control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.igfbp3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp3.ser")),  
    value)
Control: contr

REML criterion at convergence: 186.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4805 -0.6452 -0.1381  0.6411  1.7739 

Random effects:
 Groups        Name        Variance  Std.Dev. 
 litter.no:sow (Intercept) 4.980e-01 7.057e-01
 sow           (Intercept) 2.263e-18 1.504e-09
 Residual                  3.266e-01 5.715e-01
Number of obs: 72, groups:  litter.no:sow, 36; sow, 15

Fixed effects:
              Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)  7.1029694  0.4448591 33.9999998  15.967   <2e-16 ***
insem.age   -0.0001703  0.0009970 33.9999998  -0.171    0.865    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.952
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.igfbp3, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age   0.01    0.01     1    34   0.029  0.865
car::Anova(mod.igfbp3,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
           Chisq Df Pr(>Chisq)
insem.age 0.0292  1     0.8644
car::Anova(mod.igfbp3,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
               F Df Df.res Pr(>F)
insem.age 0.0281  1 27.386 0.8681

Model diagnostics

performance::check_model(mod.igfbp3)

Plot

plot.igfbp3 <- dat.l  |> 
  dplyr::filter(parameter == "igfbp3.ser") |> 
  droplevels() |> 
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 6500), 
                     breaks = seq(0, 6500, 2000),
                     minor_breaks = seq(0, 6500, by = 200) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "IGFBP3 (serum) [ng/ml]",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.igfbp3

IGFBP2/IGFBP3 (serum)

Model

mod.igfbp23 <- lmerTest::lmer(log(value) ~ 
                               insem.age +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "igfbp23.ser") |> 
                               droplevels() |>
                               drop_na(value),
                             REML    = TRUE,
                             control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.igfbp23)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp23.ser")),  
    value)
Control: contr

REML criterion at convergence: 174.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0444 -0.5111  0.1475  0.5307  1.6137 

Random effects:
 Groups        Name        Variance  Std.Dev. 
 litter.no:sow (Intercept) 4.240e-01 6.512e-01
 sow           (Intercept) 7.539e-17 8.682e-09
 Residual                  2.736e-01 5.230e-01
Number of obs: 72, groups:  litter.no:sow, 36; sow, 15

Fixed effects:
              Estimate Std. Error         df t value Pr(>|t|)   
(Intercept) -1.2353415  0.4096735 34.0000004  -3.015  0.00483 **
insem.age    0.0010469  0.0009182 34.0000004   1.140  0.26215   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.952
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.igfbp23, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age  0.356   0.356     1    34     1.3  0.262
car::Anova(mod.igfbp23,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
           Chisq Df Pr(>Chisq)
insem.age 1.3001  1     0.2542
car::Anova(mod.igfbp2_3,
           test.statistic = "F",
           type = 2)
Error: Objekt 'mod.igfbp2_3' nicht gefunden

Model diagnostics

performance::check_model(mod.igfbp23)

Plot

plot.igfbp23 <- dat.l  |> 
  dplyr::filter(parameter == "igfbp23.ser") |> 
  droplevels() |>
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 2.5), 
                     breaks = seq(0, 2.5, 0.5),
                     minor_breaks = seq(0, 2.5, by = 0.1) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "Molar ratio IGFBP2/IGFBP3 (serum)",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.igfbp23

Proteolytic activity (serum)

Model

mod.prot <- lmerTest::lmer(log(value+1) ~ 
                               insem.age +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                              dplyr::filter(parameter == "proteolysis") |> 
                              droplevels() |> 
                              drop_na(value),
                             REML    = TRUE,
                             control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.prot)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value + 1) ~ insem.age + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "proteolysis")),  
    value)
Control: contr

REML criterion at convergence: 132.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5185 -1.2360  0.2984  0.8499  1.3739 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 0.000    0.000   
 sow           (Intercept) 0.000    0.000   
 Residual                  1.339    1.157   
Number of obs: 39, groups:  litter.no:sow, 20; sow, 15

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)  
(Intercept) 1.205e+00  5.616e-01 3.700e+01   2.145   0.0386 *
insem.age   8.813e-04  1.406e-03 3.700e+01   0.627   0.5345  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.944
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.prot, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value + 1) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age  0.526   0.526     1    37   0.393  0.535
car::Anova(mod.prot,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value + 1)
           Chisq Df Pr(>Chisq)
insem.age 0.3931  1     0.5307
car::Anova(mod.prot,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value + 1)
               F Df Df.res Pr(>F)
insem.age 0.3027  1 15.645   0.59

Model diagnostics

performance::check_model(mod.prot)

Plot

plot.prot <- dat.l  |> 
  dplyr::filter(parameter == "proteolysis") |> 
  droplevels() |> 
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 21), 
                     breaks = seq(0, 21, 5),
                     minor_breaks = seq(0, 21, by = 1) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "Proteolytic activity (serum) [%]",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.prot

Stanniocalcin, STC1 (serum)

Model

mod.stc1 <- lmerTest::lmer(log(value) ~ 
                               insem.age +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "stc1") |> 
                               droplevels() |> 
                               drop_na(value),
                             REML    = TRUE,
                             control = contr)
summary(mod.stc1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "stc1")),  
    value)
Control: contr

REML criterion at convergence: 105.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.59654 -0.47741 -0.00781  0.45021  1.56136 

Random effects:
 Groups        Name        Variance Std.Dev.
 litter.no:sow (Intercept) 2.0946   1.447   
 sow           (Intercept) 1.6464   1.283   
 Residual                  0.0882   0.297   
Number of obs: 38, groups:  litter.no:sow, 19; sow, 14

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)  
(Intercept) 2.264645   1.796520 5.694156   1.261   0.2567  
insem.age   0.012553   0.005163 5.230989   2.431   0.0571 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.963
round(drop1(mod.stc1, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF DenDF F value Pr(>F)  
insem.age  0.521   0.521     1 5.231   5.911  0.057 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.stc1,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
           Chisq Df Pr(>Chisq)  
insem.age 5.9111  1    0.01505 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(mod.stc1,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
               F Df Df.res  Pr(>F)  
insem.age 4.8382  1 7.3349 0.06203 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model diagnostics

performance::check_model(mod.stc1)

Plot

plot.stc1 <- dat.l  |> 
  dplyr::filter(parameter == "stc1") |> 
  droplevels() |> 
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 10000), 
                     breaks = seq(0, 10000, 2000),
                     minor_breaks = seq(0, 10000, by = 400) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "STC1 (serum) [ng/ml]",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.stc1

Calcium (serum)

Model

mod.calc <- lmerTest::lmer(log(value) ~ 
                               insem.age +
                               # random intercept for sows and for each litter within sow
                               (1 | sow/litter.no),
                             data    = dat.l |> 
                               dplyr::filter(parameter == "calc") |> 
                               dplyr::filter(time != "30dpc") |> 
                               droplevels() |> 
                               drop_na(value),
                             REML    = TRUE,
                             control = contr)
boundary (singular) fit: see help('isSingular')
summary(mod.calc)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
   Data: drop_na(droplevels(dplyr::filter(dplyr::filter(dat.l, parameter ==  
    "calc"), time != "30dpc")), value)
Control: contr

REML criterion at convergence: 11.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.9541 -0.2379  0.0437  0.3503  1.7735 

Random effects:
 Groups        Name        Variance  Std.Dev. 
 litter.no:sow (Intercept) 4.477e-19 6.691e-10
 sow           (Intercept) 3.601e-02 1.898e-01
 Residual                  3.934e-02 1.983e-01
Number of obs: 80, groups:  litter.no:sow, 40; sow, 15

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept) 2.388e+00  9.118e-02 6.227e+01  26.190   <2e-16 ***
insem.age   2.126e-05  1.720e-04 6.730e+01   0.124    0.902    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr)
insem.age -0.804
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
round(drop1(mod.calc, test = 'Chisq'), 3)
Single term deletions using Satterthwaite's method:

Model:
log(value) ~ insem.age + (1 | sow/litter.no)
          Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
insem.age  0.001   0.001     1 67.303   0.015  0.902
car::Anova(mod.calc,
           test.statistic = "Chisq",
           type = 2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(value)
           Chisq Df Pr(>Chisq)
insem.age 0.0153  1     0.9016
car::Anova(mod.calc,
           test.statistic = "F",
           type = 2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: log(value)
               F Df Df.res Pr(>F)
insem.age 0.0151  1  25.64 0.9032

Model diagnostics

performance::check_model(mod.calc)

Plot

plot.calc <- dat.l  |> 
  dplyr::filter(parameter == "calc") |> 
  dplyr::filter(time != "30dpc") |> 
  droplevels() |> 
  drop_na(value) |> 
  # make plot
  mutate(jit = jitter(as.numeric(insem.age), 10)) |>  
  ggplot(aes(y   = value)) +
  geom_point(aes(x   = jit, 
                 col = husbandry, 
                 shape = time),
             size = 3) +
  scale_colour_manual(labels = c("Conventional",
                               "Ecological"),
                    values = c("#FFA040", "#008000")) +
  scale_y_continuous(lim = c(0, 21), 
                     breaks = seq(0, 20, 5),
                     minor_breaks = seq(0, 20, by = 2) ) +
  labs(x = "Insemination age of the sows [days]",
       y = "Calcium (serum) [mg/dl]",
       col = "Husbandry",
       shape = "Time") +
  my_theme +
  theme(legend.position = "top") +
  guides(y = guide_axis(minor.ticks = TRUE),
         col = guide_legend(override.aes = list(size  = 5)))
plot.calc

Combined plots: Figure

# Combine plots with a designated area for the legend
combined <- (plot.cort1   +
             plot.cort2   +
             plot.bioact  +
             plot.igf1    +
             plot.igf2    +
             plot.igfbp2  +
             plot.igfbp3  +
             plot.igfbp23 +
             plot.prot    +
             plot.stc1    +
             plot.calc    +
             guide_area()) + 
  plot_layout(guides = "collect") +
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(face = "bold", size = 20),
        legend.position = "right",
        legend.direction = "vertical")
combined

png("./plots/figure-age.png",
     width = 600, height = 500, units = "mm",
     pointsize = 10, res = 600)

combined

dev.off()
png 
  2 

How to cite R

“All analyses were performed using R Statistical Software (version 4.4.2; R Core Team 2024)”.

Reference: R Core Team (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

citation()
To cite R in publications use:

  R Core Team (2024). _R: A Language and Environment for Statistical
  Computing_. R Foundation for Statistical Computing, Vienna, Austria.
  <https://www.R-project.org/>.

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Manual{,
    title = {R: A Language and Environment for Statistical Computing},
    author = {{R Core Team}},
    organization = {R Foundation for Statistical Computing},
    address = {Vienna, Austria},
    year = {2024},
    url = {https://www.R-project.org/},
  }

We have invested a lot of time and effort in creating R, please cite it
when using it for data analysis. See also 'citation("pkgname")' for
citing R packages.
version$version.string
[1] "R version 4.4.2 (2024-10-31 ucrt)"
citation("performance")
Um Paket 'performance' in Publikationen zu zitieren, nutzen Sie bitte:

  Lüdecke et al., (2021). performance: An R Package for Assessment,
  Comparison and Testing of Statistical Models. Journal of Open Source
  Software, 6(60), 3139. https://doi.org/10.21105/joss.03139

Ein BibTeX-Eintrag für LaTeX-Benutzer ist

  @Article{,
    title = {{performance}: An {R} Package for Assessment, Comparison and Testing of Statistical Models},
    author = {Daniel Lüdecke and Mattan S. Ben-Shachar and Indrajeet Patil and Philip Waggoner and Dominique Makowski},
    year = {2021},
    journal = {Journal of Open Source Software},
    volume = {6},
    number = {60},
    pages = {3139},
    doi = {10.21105/joss.03139},
  }

Session Info

sessionInfo()
R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)

Matrix products: default


locale:
[1] LC_COLLATE=German_Germany.utf8  LC_CTYPE=German_Germany.utf8   
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.utf8    

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] performance_0.13.0 car_3.1-3          carData_3.0-5      lmerTest_3.1-3    
 [5] lme4_1.1-36        Matrix_1.7-1       patchwork_1.3.0    kableExtra_1.4.0  
 [9] lubridate_1.9.4    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4       
[13] purrr_1.0.2        readr_2.1.5        tidyr_1.3.1        tibble_3.2.1      
[17] ggplot2_3.5.1      tidyverse_2.0.0   

loaded via a namespace (and not attached):
 [1] gtable_0.3.6        bayestestR_0.15.1   xfun_0.50          
 [4] htmlwidgets_1.6.4   ggrepel_0.9.6       insight_1.0.1      
 [7] lattice_0.22-6      tzdb_0.4.0          numDeriv_2016.8-1.1
[10] vctrs_0.6.5         tools_4.4.2         Rdpack_2.6.2       
[13] generics_0.1.3      datawizard_1.0.0    parallel_4.4.2     
[16] pbkrtest_0.5.3      pkgconfig_2.0.3     lifecycle_1.0.4    
[19] compiler_4.4.2      farver_2.1.2        munsell_0.5.1      
[22] htmltools_0.5.8.1   yaml_2.3.10         Formula_1.2-5      
[25] pillar_1.10.1       nloptr_2.1.1        MASS_7.3-61        
[28] reformulas_0.4.0    boot_1.3-31         abind_1.4-8        
[31] nlme_3.1-166        tidyselect_1.2.1    digest_0.6.37      
[34] stringi_1.8.4       labeling_0.4.3      splines_4.4.2      
[37] fastmap_1.2.0       grid_4.4.2          colorspace_2.1-1   
[40] cli_3.6.3           magrittr_2.0.3      broom_1.0.7        
[43] withr_3.0.2         backports_1.5.0     scales_1.3.0       
[46] timechange_0.3.0    rmarkdown_2.29      hms_1.1.3          
[49] evaluate_1.0.3      knitr_1.49          rbibutils_2.3      
[52] viridisLite_0.4.2   mgcv_1.9-1          rlang_1.1.4        
[55] Rcpp_1.0.13-1       glue_1.8.0          xml2_1.3.6         
[58] see_0.9.0           svglite_2.1.3       rstudioapi_0.17.1  
[61] minqa_1.2.8         jsonlite_1.8.9      R6_2.6.1           
[64] systemfonts_1.1.0